% code to compute MT trfn from raw ts files 
% Support RTS files created by GMS05 
% usage -> mtprocn 
% uses a suit of sub-functions.

% references-----------
% 1) Sims, 1971 (ls impedance) 2) Dekker, (predicted coherency) 1984 
% 3) Gamble (rr) , 1979 (rr)	4)ProcMT manual (rts file structue, biv error)
% 5) Manoj & Nagarajan, 2003 (ANN selection)
% 6) Chave et al 1989 (Jack Knife) 7) Ritter et al 1998 (Robust processing)
% developement history-----
% 8/2000 - compute cal. table according to gain settings etc;
% 11.7.01 -computes cross and auto spectra & average it 
% 12.7 - MT app rest & phse are now computed
% 8/02 - e-pred coherency, and Biv error are calculated
% 10.2002 - rr processing introduced
% 12.2002 - rr error is being introduced
% 04.03 -- Jack knife impedance & variance
% 05.03.03 ann selection for ts segments
% 10.03.03 robspm Robust averaging for SPMs
% 15.03.03 normalize SPMs wrt power in horizontal magnetic fileds
% Robust band average and 10.6.3
% Robust section average using Jack knife 10.6.3
% Feb 22, 2011 mtprocn_h.m 

% ----known bugs & fixings
% phase distortion for 2nd and 3rd band- 05.07.2002
% OK for band1 & 4 on 07.07.2002
% 10.02 fixed phase distortion (problem was with complex transpose)
% 11.02 band2 & band1 is recorded as slices of 1024 length
%       this is not taken care of (nfft & overalp) 
% 11.02 computed error for rr is still local one
% 12.02 band1 phase distortion - problem with bd1 sys trfn ?

% ------future projects---
% ann despiker			depiker attempted bu with partial success 31.5.3
% phase rotation for time synchronization problems (?)  % Jones 1983 gave an equation for this but not yet implimented
% write edi to be completed			% not yet
% ats plugg-in -> for ADU06 ts		% not yet implimented
% latest date 02.07.2003


%create TS.matrix

% load E:\projects\magtrans\EMSLAB\EMSL13;
% load E:\projects\magtrans\EMSLAB\EMSL04;
% load E:\projects\magtrans\EMSLAB\EMSL01;
% [c,ia,ib] = intersect(emsl04.fday,emsl13.fday);
% 
% TS.matrix = zeros([length(c),5]);
% % Exl = Hxl*Zxx + Hyl*Zxy
% % Eyl = Hxl*Zyx + Hyl*Zyy
% % 
% % TS.matrix(:,1:2) = emsl01.matrix(:, 4:5); %Copy Ex and Ey to Ex and Ey slots
% % TS.matrix(:,3:5) = emsl01.matrix(:,1:3); %Copy Hx and Hy  and Hz to respective slots
% % 
% % 
% % Hxl = Hxr*Txx + Hyr*Txy;
% % Hyl = Hxr*Tyx + Hyr*Tyy;
% % 
% TS.matrix(:,1:2) = emsl04.matrix(ia,1:2); %Copy Hx and Hy to Ex and Ey slots
% TS.matrix(:,3:5) = emsl13.matrix(ib,1:3); %Copy Hx, Hy Hz to the respective slots
% 
% TS.head.nchannel = 5;
% 
% Start remote reference

% [c1,ia1,ib1] = intersect(c,emsl01.fday); %find out fday common to the three stations
% 
% TS.matrix(:,1:2) = emsl04.matrix(ia,1:2); %Copy Hx and Hy to Ex and Ey slots
% TS.matrix(:,3:5) = emsl13.matrix(ib,1:3); %Copy Hx, Hy Hz to the respective slots
% TS.matrix(:,6:7) = emsl01.matrix(ib1,1:2); %Copy Hx, Hy  to Hxr and Hyr for RR processing
% TS.head.nchannel = 7;
% ProcDef.ref =1;

%End remote reference

% Hzl = Hxr*Txx + Hyr*Txy;
% Hzr = Hxr*Tyx + Hyr*Tyy;

% TS.matrix(:,1) = emsl04.matrix(ia,3); %Copy Hz Ex slot
% TS.matrix(:,2) = emsl13.matrix(ib,3); %Copy Hz Ex slot
% TS.matrix(:,3:5) = emsl13.matrix(ib,1:3); %Copy Hx, Hy Hz to the respective slots

% Hzl = Hxr*Txx + Hyr*Txy;
% Hxl = Hxr*Tyx + Hyr*Tyy;

% TS.matrix(:,1) = emsl04.matrix(ia,3); %Copy Hz to Ex slot
% TS.matrix(:,2) = emsl04.matrix(ia,1); %Copy Hx to Ey slot
% TS.matrix(:,3:5) = emsl13.matrix(ib,[3 1 2]); %Copy Hx, Hy Hz to the respective slots

%start tipper 
% load E:\projects\magtrans\EMSLAB\EMSL01;
% TS.head.nchannel = 5;
% 
% %Hzl = Hxl*Txx + Hyl*Txy (i am interested in it)
% %Exl = Hxl*Tyx + Hyl*Tyy (some dummy)
% 
% TS.matrix = emsl01.matrix(:,[3 4 1 2 5]); %Copy Hz to Ex location, Ex to Hy locatio etc
% %end tipper

% E N D  E M S L A B

% S T A R T  I P M  A N D  P P T

%start tipper 
% ipm =load('E:\projects\magtrans\IPM_PPT\ipmnov2010_XYZ.txt'); %Quiet day XYZ
% ppt =load('E:\projects\magtrans\IPM_PPT\pptnov2010_XYZ.txt');
% TS.head.nchannel = 5;
% 
% for i = 1:3,
%     L = ipm(:,i) == 99999;
%     ipm(L,i) = median(ipm(:,i));
%     ipm(:,i) = ipm(:,i) - mean(ipm(:,i));
%         L = ppt(:,i) == 99999;
%     ppt(L,i) = median(ppt(:,i));
%     ppt(:,i) = ppt(:,i) - mean(ppt(:,i));
%end;

%load  E:\projects\magtrans\IPM_PPT\ipmpptdec2010_XYZ_mean_removed; %Quiet
%day
% 
% load E:\projects\magtrans\IPM_PPT\ipmpptnov2010_XYZ_mean_removed ipm ppt fday
% %Hzl = Hxl*Txx + Hyl*Txy (i am interested in it)
% %Exl = Hxl*Tyx + Hyl*Tyy (some dummy)
% 
% TS.fday = fday;
% TS.matrix = [ipm(:,[1 2]) ppt(:,[1 2 3])]; %Copy Hz to Ex location, Ex to Hy locatio etc
% TS.head.nchannel = 5;
% 
%end tipper

% ACE JULIA TF

%load E:\projects\ace_tensor\Electric_Field_Matrix_AllDays_128 IEF_EY_MAT EEF_MAT IEF_EZ_MAT;% The other variables are IEF_EZ_SEG EEF_SEG IEF_EY_SEG TIME_SEG;
%load /home/mnair/projects/longp/Electric_Field_Matrix_AllDays_64_JULIA_only IEF_EY_MAT EEF_MAT IEF_EZ_MAT;
load /home/mnair/projects/longp/ace_champ_data_for_mtprocn_h ACE_SEG ACE_SEG_EZ TIME_SEG CHAMP_SEG F_long Cxy_long Txy_long min_time_length len;
%JULIA DATA processing

%EEF_MAT = EEF_MAT.*24.366*1e-3; %mV/m


% EEF_MAT(:,65:end) = [];
% IEF_EY_MAT(:,65:end) = [];
% IEF_EZ_MAT(:,65:end) = [];
% 

% std_eef = std(reshape(EEF_MAT',[1,numel(EEF_MAT)]));
% length_eef = length(reshape(EEF_MAT',[1,numel(EEF_MAT)]));
% std_ief = std(reshape(IEF_EZ_MAT',[1,numel(IEF_EZ_MAT)]));
% 
% TS.matrix =  [  reshape(EEF_MAT',[1,numel(EEF_MAT)]) ; ...
%                 std_eef.*randn(1,length_eef) ; ... % Random values with norm distr and std of eef
%                 reshape(IEF_EY_MAT',[1,numel(IEF_EY_MAT)]) ; ...
%                 reshape(IEF_EZ_MAT',[1,numel(IEF_EZ_MAT)]); ...
%                 std_ief.*randn(1,length_eef)]';

% CHAMP DATA processing

std_eef = std(CHAMP_SEG);
length_eef = length(CHAMP_SEG);
std_ief = std(ACE_SEG);

TS.matrix =  [  CHAMP_SEG * 1e3 ; ...
                std_eef.*randn(1,length_eef) ; ... % Random values with norm distr and std of eef
                ACE_SEG_EZ ; ... % IEF Ez
                ACE_SEG; ... %IEF Ey
                std_ief.*randn(1,length_eef)]';
% TS.matrix =  [  std_eef.*randn(1,length_eef)*1e3 ; ...
%                 std_eef.*randn(1,length_eef)*1e3 ; ... % Random values with norm distr and std of eef
%                 std_ief.*randn(1,length_eef) ;  ...
%                 std_ief.*randn(1,length_eef) ; ...
%                 std_ief.*randn(1,length_eef)]';
% 

TS.head.nchannel = 5;            
% 
%

%
ProcDef = SetDefault_h;
ProcDef.len = floor(length(TS.matrix)/ProcDef.block);
ProcDef.SelStacks = ones([1 ProcDef.len]);  
[SPMatrix] = GetSpectra_h(TS,ProcDef);


temp = size(SPMatrix); % This is required, as number of overlapped ts sections may not be just 2*len-1
ProcDef.len = temp(2);

ProcDef.Weights = ones([ProcDef.nfrq ProcDef.len]); % Weights are modified in robust averaging routine(s)

SPMatrix = ScaleSPM(SPMatrix,ProcDef);

%disp('NormSPM-> Normalize the SPMs wrt power in horz. mag fields');
%SPMatrix = NormSPM(SPMatrix,ProcDef); % Normalize the SPMs wrt power in horz. mag fields



% compute and plot tranfer functions

 % Robust Mt
%    [tf2,ProcDef,SPM2] = robspm(SPMatrix,ProcDef);
%    subplot(211);
%    hold on;
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,1)),'k*');
% 
%    subplot(212);
%    hold on;
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,3)),'k*');

  % Coherency Thrashold
   [SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,2);
   [tf2]=tf_dn_h(SPM2);
  
open('/home/mnair/projects/longp/Figures/julia_with_error_bar.fig');
hold on;
errorbar(log10((1./ProcDef.TLFreq1)./(3600)), log10(abs(tf2(:,1))),tf2(:,5),tf2(:,5),'r');
   %% tf mag dB
%    [SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,2);
%    [tf2]=tf_dn_h(SPM2);
%    


axes('Parent',fig1,'XTick',[0.1 1 10],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);
box('on');
hold('all');
xlabel('period in hours');
ylabel('Magnitude (dB)');
hold on;
plc = 'r';
Xmag = abs(tf2(:,3));          % Spectral magnitude
Xdb = 20*log10(Xmag);   % Spectral magnitude in dB

% XdbMax = max(Xdb);      % Peak dB magnitude
% Xdbn = Xdb - XdbMax;    % Normalize to 0dB peak
% 
% dBmin = -100;           % Don't show anything lower than this
% Xdbp = max(Xdbn,dBmin); % Normalized, clipped, dB mag spec

semilogx((1./(3600*ProcDef.TLFreq1)),Xdb,plc,'LineWidth',2);
%axis([0.1 10 -55 -15 ]);

Xmag = abs(tf2(:,1));          % Spectral magnitude
Xdb = 20*log10(Xmag);   % Spectral magnitude in dB

% XdbMax = max(Xdb);      % Peak dB magnitude
% Xdbn = Xdb - XdbMax;    % Normalize to 0dB peak
% 
% dBmin = -100;           % Don't show anything lower than this
% Xdbp = max(Xdbn,dBmin); % Normalized, clipped, dB mag spec

semilogx((1./(3600*ProcDef.TLFreq1)),Xdb,'k','LineWidth',2);

grid on;

%legend('IEF E_y - EEF','(IEF E_z - EEF)^2')
    
%% Phase 
% figure1=figure;
% axes('Parent',figure1,'XTick',[0.1 1 10],...
%     'XScale','log',...
%     'XMinorTick','on');
% set(gca,'FontSize',16);
% box('on');
% hold('all');
% xlabel('period in hours');
% 
% phase_data = unwrap(angle(tf2(:,3)));
% phase_data(2) = phase_data(2)-0.15;
% % phase_data(3) = phase_data(3)-0.1641;
% % phase_data(2) = phase_data(2)-0.4784;
% 
% 
% semilogx((1./(3600*ProcDef.TLFreq1)),phase_data*180/pi,'r','LineWidth',2)
% 
% %save C:\Manoj\projects\ace\20080303\test Cxx Pxy tf phase N_data season ap_lower_limit
% axis([0.1 10, -50 ,100]);
% 
% a=axis;
% yloc=a(3)-((a(4)-a(3))*0.08);
% 
% set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10]);
% set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10']);
% text(6/60,yloc,'|<-','FontSize',16);
% text(8,yloc,'->|','FontSize',16);
% text(40/60,yloc,'-> | <-','FontSize',16);
% text(10/60,yloc,'Period in minutes','FontSize',16);
% text(2,yloc,'Period in hours','FontSize',16);
% 
%    
%    disp('Warning Errors are not computed\n');
%  % Preferential Coherency Stacking
% 	[SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,1);
%    [tf2]=tf_dn_h(SPM2);
%       subplot(221);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,1)),'r*');
%    hold on;
%    subplot(222);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,2)),'r*');
%    hold on;
%    subplot(223);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,3)),'r*');
%    hold on;
%    subplot(224);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,4)),'r*');
%    hold on;
% 
%    disp('Warning Errors are not computed\n');
%  % Jack knife impedance & error 
%    [tf2,dof]=jknf(SPMatrix,ProcDef);
%       subplot(221);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,1)),'c*');
%    hold on;
%    subplot(222);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,2)),'c*');
%    hold on;
%    subplot(223);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,3)),'c*');
%    hold on;
%    subplot(224);
%    loglog(1./ProcDef.TLFreq1,abs(tf2(:,4)),'c*');
%    hold on;
%    legend('RB','CT','CS','JK','Location','SouthEast');



   

%% Predicted time series 
% Hxl = Hxr*Txx + Hyr*Txy;
% Hyl = Hxr*Tyx + Hyr*Tyy;
%tf(:,1) = Txy, tf(:,2) = Tyx, tf(:,3) = Txx, tf(:4) = Tyy

% [TIndex,len] = GetSelArray(ProcDef);
% frq = log10(fftfrq(ProcDef.block,1/(5*60)));
% [a,ia,ib] = unique(1./ProcDef.TLFreq1);
% tf3 = tf2(ia,:);
% log_tlfreq = log10(ProcDef.TLFreq1(ia));
% 
% 
% for k = 1:4,
%     
%    [s,p] = csaps(log_tlfreq,real(tf3(:,k)));
%    sn = fnxtr(s);%This methods gives superior extrapolation !
%    tf_fft_real = ppual(sn,frq);
%    [s,p] = csaps(log_tlfreq,imag(tf3(:,k)));
%    sn = fnxtr(s);
%    tf_fft_imag = ppual(sn,frq);
%    tf_fft(:,k) = complex(tf_fft_real,tf_fft_imag);
% 
% end;
% 
% for i = 1:len,
%    
% Hxl_ts = CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),1),0);
% Hxr_ts = CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),3),0);
% Hyr_ts = CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),4),0);
% 
% 
% Hxr_spec = fft(CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),3),0));
% Hyr_spec = fft(CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),4),0));
% 
% %resample the tranfer function to that of fft
% %convert the period to log
% 
% 
% 
% 
% %Hxr * Txx
% HxrTxx  = Hxr_spec(2:end).*[tf_fft(:,3).' fliplr(tf_fft(2:end,3))']';
% HyrTxy  = Hyr_spec(2:end).*[tf_fft(:,1).' fliplr(tf_fft(2:end,1))']';
% 
% Hxl_pred_ts = real(ifft([0;HxrTxx+HyrTxy])) * 2;
% Hxl_pred_ts_from_ief_ey = real(ifft([0;HxrTxx])) * 2;
% 
% 
% subplot(211)
% plot(Hxl_ts);
% hold on;
% plot(Hxl_pred_ts,'r');
% plot(Hxl_pred_ts_from_ief_ey,'k');
% 
% subplot(212);
% plot(Hxr_ts);
% hold on;
% plot(Hyr_ts,'r');
% pause;
% clf;
% end; 	
% 
% 
% 
